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Abstract 

Two simple period determination schemes are discussed. They are 
well suited to problems involving non-sinusoidal periodic phenomena 
sampled at a few irregularly spaced points. Statistical properties 
are discussed. The techniques are applied to the double-mode Cepheids 
BK Cen and TU CAS as test cases. 

Introduction 

Fourier techniques offer the "best possible" methods of spectral 
analysis in a number of respects/ but they are not well suited to 
cases with strongly non-sinusoidal variation, or cases in which only a 
few scattered data points are available. Unfortunately, this is 
exactly the situation in many physical problems, including variable 
stars. I will discuss two classical techniques that do treat these 
cases efficiently in cases involving several discrete frequencies. 

When departing from Fourier methods, invariably some type of 
sacrifice must be made. In this case "ghosts" appear at subharmonic 
frequencies of a strong spectral line. An important point in the 
analysis is to show that the unwanted subharmonics can be minimized 
and can be distinguished from the real frequency. Once understood , 
these "ghosts" can be used to good advantage as indications of the 
presence or absence of high frequency components. 

A more complete description of these techniques will appear in 
Stellingwerf (1978). 


I. THE PPM METHOD 


A discrete set of observations can be represented by two vectors 
x, and the observation times t, where the i^h observation is given by 
** ~ 2 
(x if ti) and there are N points in all (i = 1,N). Let a be the 

variance of x, given by 



E^-x ) 2 

N-l 


( 1 ) 
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where x is the mean, x = EXj/n. For any subset of x jL 's we define the 
sample variance, s 2 , exactly as in equation (1) . Suppose we have 
chosen M distinct samples, having variances s.. 2 (j = 1,M) and contain- 

The overall variance for all the samples is then 


ing n^ data points, 
given by 


2 (n .-1) s . 2 


In^-M 


( 2 ) 


as a consequence of equation (1) . 


We wish to minimize the variance of the data with respect to the 
mean light curve. Let n be a trial period, compute a phase vector 
S>£ *= tj/n-lt.j/n] ; brackets indicate integer part. Equivalently, 

4> = t mod (II) . We now pick M samples from x using the criterion that 

all the members of sample j have similar Usually the full phase 

interval (0,1) is divided into fixed bins, but the samples may be 
chosen m any way that satisfies the criterion. All points need not 
be picked, or, alternatively, a point can belong to many samples. The 
variance of these samples gives a measure of the scatter around the 
mean light curve defined by the means of the in each sample, con- 


sidered as a function of 4>. 


We define the statistic 
0 


s’ 


(3) 


where 


2 . . o 

s is given by equation (2) and a is given by equation (1) . If 

IT is not a true period, then s 2 £ o 2 and 0 £ 1, whereas if JI is a 
correct period 0 will reach a local minimum compared with neighboring 
periods, hopefully near zero. 

Since this technique seeks to minimize the dispersion of the data 
at constant phase, we will refer to it as "phase dispersion minimiza- 
°f p pM for short. Mathematically, this is a least-squares 
fitting technique, but rather than fitting to a given curve (such as a 
Fourier component) , the fit is relative to the mean curve as defined 
by the means of each bin. We simultaneously obtain the best least- 
squares light curve and the best period. The PDM technique is thus a 
Founergram" method [as discussed by Faulkner (1977)] of infinite 
order, since all harmonics are included in the fitted function. The 

? e , r f es technique, a least squares fit to a truncated series 
with variable amplitudes and phase, often requires additional con- 

1976) ntS Snd rathor high orders for nonsinusoidal variations (Lucy, 


Although the individual samples may be chosen in many ways, it is 
convenient to define a standard bin structure. We divide the unit 
interval into bins of length 1/N b , and take N Q "covers" of bins, 

each cover offset in phase by l/(N fa N c ) from the previous cover, using 
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periodic conditions on the unit interval to obtain a uniform covering. 
We thus obtain M = Nj j N c bins, each of length 1/N b , and whose midpoints 

are uniformly spaced along the unit interval at distance of 1/Oyy . 

Clearly, each data point will fall in exactly N c bins. Denote a given 

bin structure by (N^, N c ) . 

To compute 0, we take the ratio of variances of the two subsets 
of X, that of the actual observations, x, and that of the bins. 0 

therefore has a probability density given by an F distribution with 
En-M and N-l degrees of freedom. It is convenient- to define F as a 

number greater than unity, so F = 0 1 . The probability,?, that a 
given value of 0 is due to random fluctuations (also called the^ 

"significance") is twice the area of the F distribution above 0 
(two-sided test). This probability approaches unity as 0 1. Thus, 
for significance P, we compute 


'(P/2,B 1£ ,H 2£ ) - V®. H lf = n-l, l. 2f - En.-H. 


P may then be obtained by reference to an F table, or using an approxi- 
mation to P(F,N lf N 2f ), see Abramowitz and- Segun (1965), I 26. 

If N is large (> 100), we may take o^(x) £ o^(X). In this case 

2 ~ 

we may use the somewhat simpler x test. 


X 2 (P/2,N f ) = N f° ' N f = 

Here P is twice the area of the x distribution below N^G. 


(5) 


II. THE WR METHOD 

An interesting related method is discussed by Whittaker and 
Robinson (1926, "WR") in which one seeks the maximum variance of the 
bin means (as opposed to the mean of the bin variances) • If is 

the variance of the bin means, define 0^ « l-s m ^/0^ for comparison 

purposes. In general 0^^ will vary between 0 and l-<l/nj>, where <n j > 

is the mean number of points per bin. At a true period s^ 2 £ ° and 

0 =0. Here we seek periods at which the amplitude of the mean curve 

WR , , 

is a maximum, which in most cases will correspond to minimum phase 

dispersion. The calculation of is easier than s (equation 2) , but 

the number of degrees of freedom is much lower, suggesting less 
sensitivity (see § III) • We will show below that in most cases this 
supposition is true, although for one range of parameters the WR 
technique may be preferrable. The WR method is discussed in detail by 
Chapman and Bartels (1940) under the title "persistence analysis". 
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A statistical criterion may also be derived for the WR technique. 

2 

r n . . 

3 


The distribution of bin means is normal if M is large with o m 2 = o^/n 


Since M is generally not large, however, a t distribution with M— 1 
degrees of freedom will actually be obtained, with 


2 2 

°m % (M-l)/(M-3) a /ny We may always select bins with n^ = N/N^, 

all j. Note that normal constant-phase bins will produce binomially 

distributed n . and increase a 2 
2 3 m 


for 


m 

than a. 


since a harmonic mean of n^ appears in 
We wish to test whether the observed s 2 is significantly bigger 


m 


m 

• 2 2 

Evidently s ra /° m follows an F distribution with M-l and N-l 


degrees of freedom. Using the definition of we have 


WR' 


F *0 (M-3) ’ N ,, . 

(M-l, N-l)- (M-l) U-9^)-. . (6) 


* HI. APPLICATIONS 

Here I will briefly describe the practical aspects of these 
methods; for details see Stellingwerf (1978). 


Typical transforms are illustrated by Figure 
computed 0-transforms for a sine wave fsin(2Hft)] 


l f which shows the 
panel A, and a 



Fig. 1: 0 statistic versus frequency for the two test 

cases described in the text. Panel A: sine-wave 

transform; Panel B: sawtooth function transform. 
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"saw-tooth" function (fractional part of ft), panel B. These functions 
have been selected because they represent the two extremes found in 
variable star light curves. In each case x consisted of 201 data 
points evenly distributed over ten periods^with f = 1, T = 10. A bin 
structure of (N b ,N c ) = (5,2) was used. The main line at f = 1 is 

virtually identical to the Fourier power spectrum line shape 


f (sin (x)/x) 2 ] . Two subharmonics are present (the third, just visible 
at f «= 1/4 is below the cutoff n = N b £ 3, and is not significantly 

different from 1). The great similarity of the two curves shown in 
Figure 1 is, of course, the strong point of the PDM technique: highly 
nonsinusoidal variations are handled just as efficiently as a sine wave. 
The presence of subharmonic response could be a disadvantage if oscil- 
lations with widely spaced frequencies are present. In practice, 
however, subharmonics can be distinguished in at least three ways: 

1) light curve shape, 2) narrow line widths, 3) reduced significance 
with increasing bin size. If initial scans of the full frequency range 
use broad bin sizes [(5,2), say], subharmonics should pose no problems. 


It is important to know the number of observations required to^ 
achieve a given significance with these methods. This quantity. ("N") 
is plotted in Figures 2 and 3 as a function of the signal-to-noise 
ratio, e = o/^oi^e* A var * et y bi n structures are shown. It is 

seen that coarse bins and multiple covers are needed for small N and e. 
Clearly, the PDM technique is superior at small N, while the WR 
technique seems preferable in noisy cases (e < 0) if N is large enough. 
Accuracy is lost if bins are too wide, and multiple covers will not 
help if many bins contain the same points. In general (5,2) structure 
works very well. To obtain mean curves for prewhitening, bins should 
be made as fine as possible. 


IV. MULTIPLE PERIODS OF BK CEN 

To illustrate the use of these methods, we briefly present an 
analysis of the light variation of the double-mode Cepheid BK Centauri. 
For this purpose we will use only the 49 photoelectric measurements 
made by C. J. van Houten in 1965 for which Leotta-Janin (1967) has 
published the AV values. This author comments that "their number is 
too small to allow a reliable determination of the beat period" and so 
uses 25 years of plate estimates as well as knowledge of other Cepheid 

period ratios to obtain Jig = 3^17389, 11^ = 2f2366. We have found that 

the photoelectric points alone unambiguously confirm these results, and 
provide information on mode interaction as well. 

The 49 measures span 137 days. Line widths will therefore be 
about 0.007 in frequency, so the frequency step size was taken to be 
.0015 for moderate resolution. Twelve nights have two observations, 
providing minimal alias discrimination. Although the measurements are 
quite accurate, we nonetheless estimate e 1.5 due to the secondary 
oscillation. Figures 2 and 3 indicate that the PDM and the WR 
techniques should be about equally good for this case, provided Nj^lO. 
We choose a (5,2) bin structure to maintain about 10 points per bin. 
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Figure 4, panel A shows the 0 transform of the raw data. The 
primary period is clearly evident (minimum e) , with minima a and b 
forming the standard subharmonic sequence found in Figure 1. Minimum 
f is narrow and showed a -two-cycle mean curve, it is therefore the 



Figure 4 s 0 transforms of the BK Cen data, plotted versus 
frequency in cycles/day. Panel A: transform of 

the raw data. Panel B: transform of the data 

after removal of the fundamental oscillation. 

Panel C: transform of the data after removal of 

the fundamental and first overtone oscillations. 
Labeled minima are identified in Table 1. 

first subharmonic of a sizeable minimum off the scale, at f % 0.7. 

This is the main alias of e, occuring at 1-fp. Comparison of the sub- 
harmonics band f indicates that e is the principal frequency. The 
identifications of the minima, together with their approximate fre- 
quencies and significances, are given in Table 1. Among the barely 
significant group, we identify as the second subharmonic of the 
second alias of e, at l+f«, while d, h and i are first overtone 
features. w “ “ “ 

We see here a distinct advantage of the subharmonic response of 
this methods any other major frequency in the range . 6<f<1.2 would 
show a subharmonic on Figure 1(A) comparable to b or f. Any minimum 
in the range 1.2<f<1.8 would show a second subharmonic comparable to 
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TABLE 1 


Identification of Features in Figure 4 

Feature 

Identification 

Frequency (d“ A ) 

Significance 

a 

v j 

.105 

0.13 

b 

V 2 

• 158 

0.00 

c 

(l-f 0 )/3 

.226 

0.39 

A 

a - fj>/2 

.276 

0.09 

ft 

*0 

• 316 

0.00 

f 

U - f 0 >/2 

.344 

0.00 

9 

a + t 0 )/3 

• 439 

0.24 

h 


.447 

0.27 

i 


• 553 

0.05 

3 

V 3 

*• 146 

0.50 

Sc 

(l-f 1 )/3 

.103 

0.55 

X 

V 2 

.224 

0.04 

n 

( l - fjJ /2 

.276 

0.08 

n 

(f 0 +fi>/ 2 

• 385 

0.33 

o 

*1 

.447 

0.00 

P 

w* 

• 553 

0.00 

n 

Cl-Jfo+fjJJ/Z 

.116 

0.60 

r 


.235 

0.27 

8 

{ VfiJ/3 

.253 

0.55 

t 

{ f 0 +£ x )/2 

.383 

0.43 

u 


• 764 

0.11 


a or g. Since all features have been identified, no high frequency 
components exist. This "look-ahead" feature can be extended by in- 
creasing N^, at the expense of further complicating the spectrum. 

Having identified the principal frequency (fg) , this component 

was then removed from the data (using linear interpolation between bin 
means) . The transform of the reduced data is shown in panel B of 
Figure 1. Here j_, 1, o are due to the first overtone, while k, m, £ 
represent the aliases. Note that each alias feature is slightly less 
significant than the corresponding "real" feature. A mode interaction 
feature, n, has also appeared. We can again conclude that no more 
significant modes exist to f = 1.8, since no subharmonics appear. 

Panel C shows the transform of the doubly- reduced data. No 
significant feature appears, but many marginal features are related to 
the mode interaction term fg+f^. An extension of the transform to 

higher frequencies showed that this mode did indeed appear at the 
P = 0.11 level. This is the only remaining significant mode. This 
particular interaction term (f Q +f^) has also been found to be excited 

in several other studies of similar objects (Fitch 1976, Henden 1976, 
Faulkner 1977a, b, Fitch et.al. 1978). In particular, this term invari- 
ably dominates the beat frequency (f^-fp) component. At present the 
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meaning of this is not clear. 


The periods and amplitudes of the three components were determined 
using a finer frequency scan, and are given in Table 2. Allowing for 
the bin size, the actual amplitudes should be about 20% larger, this is 

TABLE 2 


Summary of Principal Components 


Mode 

Frequency ^ * 

Period ^ 

Amplitude 

Adjusted 

Amplitude 

Phase* 

f 0 

0.3153 

3.172 

0?S44 

0?C48 

0.40 * 

1 1 

0.4473 

2.236 

0?222 

0?204 

0.56 

f D +f l 

0.7602 

1.302 

O.'lSl 

0?180 

0.16 


•Phase - ttjnax-tiJ/V *i “ «>2, 438. 813.48 

shown in the "corrected amplitudes". The sum of these estimates is 
1.1, in agreement with the range of the photoelectric measures. The 
standard deviation of the residuals (triply-reduced data) was 0^048 . 

The frequency labeled f 0 +f x in Table 2 is listed as derived and 

the sum of the principal components. The discrepancy may 
not be real, since it amounts to about 1/T — suggesting a side lobe 
problem. 
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Discussion 


A. Cox: What is the a of the data after you've taken the periods out? 

Stellingwerf : The residual a following the triple reduction was a « 0.05 

magnitudes. 
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